
#|____________________________________________________________________
 |
 |         RMDS.LSP - Replicated Multidimensional Scaling
 |              Copyright (c) 1997, Forrest W. Young
 |                        SMACOF ALGORITHM
 |____________________________________________________________________
 |#



(defmeth rmds-model-proto :matrices-to-list (mats) 
  (let* ((m (length mats))
         (n (send self :nstim))
         (out-list (repeat nil (/ (* m n (- n 1)) 2)))
         (k -1)
         )
    (dotimes (matk m)
             (dotimes (i (- n 1))
                      (dolist (j (iseq (+ i 1) (- n 1)))
                              (setf k (+ k 1))
                              (setf (select out-list k) 
                                    (aref (select mats matk) i j)))))
    out-list))

(defmeth rmds-model-proto :torgmds (O R)
"Args: O R
Computes Torgerson's Classical MDS of observed data in O, a symmetric (NxN) matrix of dissimilarities. Returns  an an (NxR) matrix of R-dimensional coordinates of points in Euclidean space. Optimizes fit of XX' to B, where B is double-centered (^O 2)."
  (let* ((osq (^ o 2))
         (n (send self :nstim))
         (means (mapcar #'mean (row-list osq)))
         (ones (repeat 1 n))
         (col-means (outer-product ones means))
         (B (* -.5 (+ (- osq col-means (transpose col-means)) (mean means))))
         (svd (sv-decomp B))
         (singvecs (first svd))
         (singvals (second svd))
         (L (diagonal (sqrt (select singvals (iseq r)))))
         (U (select (first svd) (iseq n) (iseq r))))
    (matmult U L)))


(defmeth rmds-model-proto :guttman-update (X O D)
"Args: X O D
Where X is an nxr matrix, O is a list of m nxn matrices, and D is an nxn matrix. Returns X+, Guttman's update of the coordinates X by computing X+ = BX, where the matrix B is based on observed data O and the current distances D amongst X."
  (let* ((n (send self :nstim))
         (m (send self :nmats))
         (D+I (+ D (identity-matrix n))) 
         (OMean (/ (apply #'+ O) m))
         (O/D (- (/ OMean D+I)))
         (sums (mapcar #'sum (row-list O/D)))
         (B (- O/D (diagonal sums))))
    (/ (matmult B X) n)))


(defmeth rmds-model-proto :r-stress (D O-List)
"Args: D O-List
D is an nxn matrix. O is a list of m nxn matrices. Computes Kruskal's replicated stress between the distances D and the observed data in O-List."
  (let* ((s2 0)
         (m (send self :nmats)))
    (dotimes (i m)
             (setf s2 (+ s2 (^ (send self :stress D (select O-List i)) 2))))
    (sqrt (/ s2 m))))

(defmeth rmds-model-proto :stress (D O)
"Args: D O
D is an nxn matrix. O is an nxn matrix. Computes Kruskal's stress between the distances D and the observed data O."
  (sqrt (/ (sum (^ (- D O) 2)) (sum (^ O 2)))))

(defmeth rmds-model-proto :zero-diagonal (mat-list)
"Args: MAT-LIST
MAT-LIST is a list of square matrices. Sets diagonal of each matrix to zero."
  (mapcar #'(lambda (X) (- X (diagonal (diagonal X)))) MAT-LIST))
  

(defmeth rmds-model-proto :norm-diss-data (mat-list)
"Args: mat-list
MAT-LIST is a list of square nxn matrices. Normalizes each matrix so that lower triangle has rms 1"
  (let ((n (send self :nstim)))
    (mapcar #'(lambda (O) (/ O (sqrt (/ (sum (* O O)) (* 2 n (- n 1))))))
                mat-list)))

(defmeth rmds-model-proto :r-lsmtransform (ord-mat-list int-mat)
"Args: ORD-MAT-LIST INT-MAT
Manages least squares monotonic transformations for replicated MDS. ORD-MAT-LIST is a list of m nxn matrices of ordinal data. INT-MAT is an nxn matrix of quantitative distances. Returns a list with four elements: 1) List of matrices with transformed data; 2) list of ordered data; 3) list of transformed data; 4) list of ordered data"
  (let* ((m (send self :nmats))
         (list-of-trf-matrices nil)
         (ordered-ord-list nil)
         (ordered-data-lists nil)
         (ordered-trf-list nil)
         (ordered-int-list nil)
         (result nil)
         )
    (dotimes (mat m)
             (setf result (send self :lsmtransform (select ord-mat-list mat) int-mat))
             (setf list-of-trf-matrices 
                   (make-matrix-list list-of-trf-matrices (select result 0)))
             (setf ordered-ord-list (append ordered-ord-list
                                            (list (select result 1))))
             (setf ordered-trf-list (append ordered-trf-list
                                            (list (select result 2))))
             (setf ordered-int-list (append ordered-int-list
                                            (list (select result 3))))
             )
    (list list-of-trf-matrices 
          ordered-ord-list
          ordered-trf-list
          ordered-int-list)))


(defmeth rmds-model-proto :lsmtransform (ord-mat int-mat)
"Args: ORD-MAT INT-MAT
Reshapes ORD-MAT (ordinal data matrix) and INT-MAT (interval distance matrix) from symmetric matrices to lists, using only lower triangular elements. Then uses LSMT to compute transform using lists. Returns a list whose elements are 1) the least squares monotonic tranformed data in symmetric matrix form; 2) the data in ordered list form; 3) the transformed data in ordered list form; and 4) the distances in ordered list form."
  (let* ((n (send self :nstim))
         (ord-list (repeat nil (/ (* n (- n 1)) 2)))
         (int-list (repeat nil (/ (* n (- n 1)) 2)))
         (trf-list nil)
         (ordered-ord-list nil)
         (ordered-trf-list nil)
         (ordered-int-list nil)
         (trf-mat (make-array (list n n) :initial-element 0))
         (ranks    nil)
         (k -1)
         )
    (dotimes (i (- n 1))
             (dolist (j (iseq (+ i 1) (- n 1)))
                     (setf k (+ k 1))
                     (setf (select ord-list k) (aref ord-mat i j))
                     (setf (select int-list k) (aref int-mat i j))))
    (setf trf-list (lsmt ord-list int-list))
    (setf trf-list (* trf-list (/ (ssq int-list) (ssq trf-list))))
    (setf k -1)
    (dotimes (i (- n 1))
             (dolist (j (iseq (+ i 1) (- n 1)))
                     (setf k (+ k 1))
                     (setf (aref trf-mat i j) (select trf-list k))
                     (setf (aref trf-mat j i) (select trf-list k))))
    (setf ranks (order ord-list))
    (setf ordered-ord-list (select ord-list ranks))
    (setf ordered-trf-list (select trf-list ranks))
    (setf ordered-int-list (select int-list ranks))
    (list trf-mat ordered-ord-list ordered-trf-list ordered-int-list)))
                           
;; =-=-=-=-=-=-=-=--=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
;; lsmt function is in regress.lsp. Should be moved to function.lsp
;; would be nice to add primary ties approach
;; =-=-=-=-=-=-=-=--=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

(defun lsmt (ordinal interval)
"Args: ORDINAL INTERVAL
Computes the least squares monotonic transformation of ORDINAL given INTERVAL.
ORDINAL  is a vector or list of numbers of an ordinal  variable.
INTERVAL is a vector or list of numbers of an interval variable.
ORDINAL and INTERVAL must have the same number of elements.
Returns a vector or list of the numbers which are (weakly) in the same order
as ORDINAL and which are a least squares fit to INTERVAL. 
The secondary approach to ties is used (ties remain tied).
Written by Forrest W. Young 10/26/94"
(when (/= (length ordinal) (length interval))
        (error "The ordinal and interval variables must be the same length."))
  (let* ((n (length ordinal))
         (rank-order (order ordinal))
         (X (select interval rank-order))
         (Y (select ordinal  rank-order))
         (block-mean nil)
         (block-size nil)
         (tie-mean nil)
         (j nil)
         (lower nil)
         ) 
;force tied data to remain tied
    (when (/= n (length (remove-duplicates ordinal)))
          (dolist (i (iseq 0 (- n 2)))
                  (when (and (= (select Y i) (select Y (1+ i))) 
                             (not lower)) 
                        (setf lower i))
                  (when (and lower 
                             (or (< (select Y i) (select Y (1+ i)))
                                 (= i (- n 2))))
                        (setf tie-mean 0)
                        (when (= i (- n 2)) (setf i (- n 1)))
                        (dolist (j (iseq lower i))
                                (setf tie-mean (+ tie-mean (select X j))))
                        (setf tie-mean (/ tie-mean (- (1+ i) lower)))
                        (dolist (j (iseq lower i))
                                (setf (select X j) tie-mean))
                        (setf lower nil))))
    (setf Y (copy-list X))
;compute the least squares monotonic transformation
    (dolist (i (iseq 1 (1- n)))
            (setf block-size 1)
            (setf block-mean (select x i))
            (loop 
             (setf j (- i block-size))
             (when (or (< j 0) (>= block-mean (select y j))) (return))
             (setf block-size (1+ block-size))
             (setf block-mean  (/ (+ (* (1- block-size) block-mean) 
                                     (select x j)) block-size)))
            (dolist (k (iseq (1+ j) i)) (setf (select Y k) block-mean)))
    (select Y (order rank-order))))

#| fwy changed 09-28-02

(defun get-data-matrices (dob)
  (let* ((n (send dob :active-nvar '(numeric)))
         (m (send dob :active-nmat '(symmetric)))
         (data (send dob :active-data-matrix))
         (data-matrix-list (list (matrix (list n n) (col data 0))))
         )
    (dotimes (i (- m 1))
             (setf data-matrix-list 
                   (make-matrix-list data-matrix-list
                                     (matrix (list n n) (col data (+ i 1))))))
    data-matrix-list))
|#

(defun get-data-matrices (dob)
  (send dob :get-active-data-matrices))


